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Abstract 


Q 


Lift coefficient 


The nonlinear indicial response method is used to 
model the unsteady aerodynamic coefficients in the low 
speed longitudinal oscillatory wind tunnel test data of 
the 0.1 scale model of the F-16XL aircraft. Exponential 
functions are used to approximate the deficiency func- 
tion in the indicial response. Using one set of oscilla- 
tory wind tunnel data and parameter identification 
method, the unknown parameters in the exponential 
functions are estimated. The genetic algorithm is used 
as a least square minimizing algorithm. The assumed 
model structures and parameter estimates are validated 
by comparing the predictions with other sets of avail- 
able oscillatory wind tunnel test data. 

Nomenclature 




a(a),a*(OL),a 2 *(OL) Indicial parameters in Model I and 
Model II 


dC 

Lift-Curve-Slope, — — , per rad 
da 

Rotary derivative in pitch, 

C„ = C L or C nr per rad 



Pitching moment coefficient 

Slope of pitching-moment-coefficient 
dC 

curve, — — , per rad 
da 


Acceleration derivative in pitch, 


dC 


3 ^, per rad 


ac 

2V 


b i (a) i b ] '(a),b 2 '(a) Indicial parameters (exponents) in 
Model I and II, 1/sec 


F 2 , F }J F a 

J 


Deficiency functions in Model I and II 
Cost function as defined in Eq. (23) 


C| ,C 2 ,C 3 


Constants in the algebraic model 
for 


k 


Reduced frequency, 


cor 

2V 


Mean aerodynamic chord (mac), m 


M 


Mach number 




Constants in the algebraic model 
for C mq 


Pitch rate, rad/sec 
Time, sec 
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Period of oscillation, sec 
Velocity, m/sec 

Distance to moment reference point (in 
terms of mac) from nose, m 
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Distance lo center of the planform area 
(in terms of mac) from nose, m 


Angle of attack, deg 



(0 Frequency of oscillatory motion in 

pitch, rad/sec 

x Dummy variable for time 

Acronym 

GA Genetic Algorithm 

Introduction 

With the advent of innovative high alpha control 
devices such as three-axes thrust vectoring and fore- 
body vortex controllers* it has become increasingly 
attractive for fighter aircraft to exploit the full potential 
offered by flight at high angles of attack (at or beyond 
stall) to achieve air superiority in close-air combat. 
Since aircraft development programs tend to rely heav- 
ily on simulations to reduce the risk and cost of flight 
test programs, accurate modeling of such motions as- 
sumes critical importance. 1 The traditional aerodynamic 
approach based on time invariant stability and control 
derivatives proves to be inadequate for predicting such 
aircraft motions because it does not account for un- 
steady aerodynamic effects that are significant at such 
flight conditions. The problem of modeling unsteady 
aerodynamic effects at high angles of attack is a formi- 
dable task because the flow field is characterized by 
extensive flow separations, vortices, vortex interactions 
and possible vortex bursts over the wing and tail sur- 
faces. In view of this, generalized unsteady aerody- 
namic models that are suitable for high angles of attack 
flight simulation are not yet available. 

NASA Langley Research Center has conducted 
low speed static and forced oscillation tests on the 0.1 
scale model of the F-16XL aircraft at high angles of 
attack in the Langley 12-Foot Low Speed Wind Tunnel. 
Using these experimental data and the indicial response 
method, 2 * 3 an attempt has been made in this paper to 
develop aerodynamic models suitable for predicting 
unsteady aerodynamic coefficients in the longitudinal 
forced oscillatory motion of the F-16XL aircraft model 
at high angles of attack. 

The indicial response method provides a funda- 
mental approach to the problem of modeling unsteady 
aerodynamic effects. One significant advantage of this 
approach is that when the indicial response to a par- 
ticular forcing mode is known, e.g., due to angle of 
attack, the cumulative response to an arbitrary forcing 
mode can be obtained using the superposition principle 
in the form of Duhamel integral. By definition, an indi- 
cial function is the response to a change in the surface 


boundary condition, say angle of attack, that is applied 
instantaneously and held constant thereafter. The indi- 
cial response is a mathematical concept and as such is 
difficult to determine by an experiment. In general, the 
indicial response can be assumed to consist of two dis- 
tinct components, the noncirculatory loading and the 
circulatory loading. The noncirculatory loading arises 
due to the sudden change in the surface boundary con- 
dition such as a sudden change in the angle of attack. It 
assumes a large peak value (theoretically infinity in 
incompressible flow) initially but subsequently decays 
rapidly. At present, results based on the linear piston 
theory are used to define the initial value of the noncir- 
culatory loading. For two-dimensional airfoils in in- 
compressible flow, Lomax 4 has presented analytical 
expressions for the early part of the decay of the 
noncirculatory loading. For arbitrary configurations, 
such information is not available. The circulatory load- 
ing is initially zero and builds up gradually to its final 
steady state value that can be obtained in a relatively 
easier fashion either by CFD methods or static wind 
tunnel tests. The difficulty lies in modeling the decay of 
noncirculatory loading and the build up of the circula- 
tory loading. Approximating the decay of the noncir- 
culatory loading and the build up of circulatory loading 
with exponential functions and using linear indicial 
response method, Beddoes, 5 and, Leishman 6 have ob- 
tained satisfactory results for the oscillatory motion of 
two-dimensional airfoils at low angles of attack in sub- 
sonic flow. In the linear indicial response method, the 
indicial response is independent of angle of attack. In 
other words, if we know the indicial response at one 
angle of attack, the same is applicable for all other an- 
gles of attack as long as aerodynamic coefficients are 
linear with angle of attack. 

In this paper, the nonlinear indicial response 
method developed by Tobak 2 3 has been used to analyze 
the low speed, longitudinal oscillatory wind tunnel data 
of the F-16XL aircraft model at high angles of attack. 
The linear indicial method is not suitable for high an- 
gles of attack because, in general, the aerodynamic co- 
efficients vary nonlinearly with angle of attack. In the 
nonlinear indicial response method, the indicial re- 
sponse depends on angle of attack. In Tobak’ s 2,3 for- 
mulation, the time varying part of the indicia! response 
is represented by the “deficiency” function. Thus, the 
task now becomes one of modeling the time history of 
the deficiency function. In this paper, exponential func- 
tions are used to approximate the time history of defi- 
ciency function. First, we use a model with two term 
exponential functions, one to approximate the decay of 
noncirculatory part and the other to approximate the 
decay of circulatory part. This two term exponential 
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representation is called model I and is physically con- 
sistent because it explicitly models the decay of both 
noncirculatory and circulatory components of the defi- 
ciency function. According to Tobak, 2 3 the initial mag- 
nitudes of the noncirculatory and circulatory compo- 
nents are fixed by the theory. Therefore* the task is to 
model the decay of noncirculatory and circulatory com- 
ponents of the decay function. 

Secondly, we use one term exponential function to 
approximate the complete deficiency function as pro- 
posed by Klein 7 ” 9 et al which is called Model IT in this 
study. This one term model has a simple structure and 
is computationally efficient because it can be expressed 
in the form of an ordinary differential equation as op- 
posed to the convolution integral form of Model L By 
comparing the predictions of Model II with those of 
Model I, it is proposed to evaluate the validity of Model 
II for predicting the longitudinal oscillatory unsteady 
aerodynamic coefficients of the F-16XL aircraft model 
at high angles of attack. 

In this paper, parameter identification is used to 
determine the unknown parameters in the exponential 
functions of Model I and Model II. Some adjustments 
have been introduced in Model I to account for uncer- 
tainties in the use of the static wind tunnel data that is 
needed in Model I and to account for deviations from 
linear piston theory due to nonlinearities at high angles 
of attack. Only one set of longitudinal oscillatory test 
data is used for parameter estimation purposes. The 
other sets of available oscillatory test data are used to 
validate the parameter estimates. The genetic algorithm 
is used to minimize the sum of squares of error between 
the estimated and measured time histories of lift or 
pitching moment coefficients. 

Analysis 

Consider the problem of estimation of unsteady lift 
and pitching moment coefficients. According to the 
nonlinear indicial function formulation of Tobak, 23 the 
time varying lift and picthing moment coefficients 
during oscillatory motion in pitch, for first order accu- 
racy in frequency, can be expressed (in Tobak’s 2 3 nota- 
tion) as 


C t (0 = C t (oo,a(f), <7 = 0) + f^jc t4 
-J F,(r - r,a(t ),<7 = 0)a(x)t/t 


( 1 ) 


cj<) = cj~Mt),q = Q)+[^X: mil 

( 2 ) 

-j F,(7 - t; a.(x),q- 0 )d( x)dx 

0 

where d = —. In this formulation, the terms like q 
dt 

and dq are ignored since they will be of second order 
in frequency. The first term on the right hand side of 
Eqs (1) is the lift coefficient in steady flow conditions 
(/ ^ oo) and at an angle of attack equal to a(t) and zero 
pitch rate. Similarly, the first term on the right hand 
side of Eq.(2) is the pitching moment coefficient under 
the steady flow conditions. Both these coefficients can 
be obtained either from static wind tunnel tests or CFD 
methods. The second terms on the right hand sides of 
Eq. (1) and (2) contain the familiar damping in pitch 
derivatives. Usually, the isolated value of or C nh( is 
not known because the small amplitude forced oscilla- 
tion wind tunnel tests that are normally used to obtain 
pitch damping derivatives give the combined value of 

< C L, + C ui) or ^ + In such lesls with r ' = a ’ il 

is not possible isolate values of C Ll( or C nhJ . In view of 
this, both C Ul and C nUj are assumed unknown. The third 
terms on the right hand sides of Eq.(l) and (2) with 
convolution integrals represent the unsteady lift and 
pitching moment coefficient respectively and contain 
the unknown deficiency functions F 2 and F 3 which de- 
pend on angle of attack and the elapsed time t - x. Ac- 
cording to Tobak, 2 3 

F 2 (a,t) = C Lu (a,t = <*>)- C La (<x, t) ( 3 ) 

F i (a,t) = C„Ja,t = c»)-C„Ja,t) ( 4 ) 

where QJoqf) and C ma (a,t) are the nonlinear, time de- 
pendent indicial responses, and C La ( aj = °o) and 
C ma ( <x,t = oo) are their steady state values. Therefore, as 
(t oo), F 2 ( a,0 = F 3 (a,r) = 0. 

Uet 


>■,(0 = C t (f)-C t (~.c«0.0)-^C I , (5) 

y 2 (0 = CJO-CJ o=,a(0,0)-^jc^ (6) 

so that 

y,(0 = -J F 2 « - x;a(x),q = 0)a(x)dx (7) 

ft 
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y 2 (t) = -J F,(/ - x;a(T ),<7 = 0)a(t)Jt (8) 

Wc observe that y ,(>) and v 2 (0 denote the unsteady 
components of time varying lift and pitching moment 
coefficient respectively. 

Model I 


like surfaces at zero angle of attack. Therefore, for ap- 
plication to complex aircraft configurations such as F- 
16XL operating at high angles of attack, it is necessary 
to derive a suitable initial value for the nonciculatory 
component. This, however was not attempted in this 
study. Instead, to account for the deviations from linear 
piston theory, an unknown empirical multiplier a 2 *(OL) is 
introduced in the second term. With these assumptions. 


We assume that deficiency functions F 2 (olj) and 
F 3 ( cu) consist of two term exponential functions as 
given by 


F 2 (a,t) = C L Ja)e- h " a) ' 

F,(a,/) = 



x rs )e 


-ftl a)/ 


(9) 

( 10 ) 


F 2 (a,t) = a l \a)C La (a)e- h ' ia) ' 

Substituting in Eqs. (7) and (8), wc obtain, 


(II) 


( 12 ) 


Here, Q„(a), and C M (a), X ‘* ) 

arc the terms fixed by the theory, 2,3 b { ( a), b x \ a), b 2 {a) 
and b 2 \a) are unknown indieial parameters which are 
assumed to be functions of angle attack,*^ is the dis- 
tance (in terms of mean aerodynamic chord, mac) from 
the nose to the moment reference point, and x cg is the 

distance (in terms of mac) from the nose to the center of 
the planform area of the vehicle. The determination of 
all unknown parameters will be discussed later. The 
first term on the right hand side of Eqs. (9) or (10) ap- 
proximates the decay of circulatory component. To 
evaluate this term, we need to know the “theoretically” 
correct value of static lift-curve-slope, Q^oc) or the 
pitching-moment-curve slope C ma ( a). However, to ac- 
count for uncertainties in determining a “theoretically” 
correct value of C^fa) or C m{X { a) using static wind tun- 
nel tests, we introduce a multiplier tf,*(oc). This un- 
known parameter is assumed to be a function of angle 
of attack (with suitable upper and lower bounds). The 
second term in Eq. (9) approximates the decay of 
noncirculatory component with an initial value equal to 

f — j. This initial value is based on the linear piston 

km) 

theory and is supposed to be independent of planform 
shape of the body. The noncirculatory component is 
assumed to be uniformly distributed over the entire 
planform surface of the body so that the resultant force 
acts at the center of the planform area giving a moment 
arm ( x reJ - x cg ) as used in Eq.(10). However, as noted 

by Tobak, 10 this result is applicable mainly for flat plate 


}j(0 


J a.Xt))^ (a(T))p M " ,l))( '"'a(TVi 
0 

-(— lja 2 ’(a(T))<’ -6(la,t))l '‘ t, d(T)</T 


(13) 


ia l ’(a(T))C, ro (a(t))e"'- (o(,>)( '' T) d(T)t/x 

0 

K M 

(14) 


Note that when the “theoretically correct” value of 
C L (a) and C m (a) are known, a) (a) = 1, and for flat 
plate like surface operating at low angles of attack, 
a 2 \a) = 1 . 

The rotary derivatives are assumed as follows: 9 

C U{ = c } , a <20 deg (15) 



20 deg < a < 70 deg 


C m =d x +dS— ) + 1 , 

m “ 1 -I57.3J \57.3; < 

0 < a < 70 deg 

where c„ c 2 , c 3 , d x , d 2i and d 3 are unknown constants. 
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Model II The model structures for C U{ and C, Ml/ as given in 

Eqs. ( 1 5) to ( i 7) for Model I were also used for Model IT. 

We approximate the complete deficiency function 
by a single exponential term 7 * 8 9 as follows: Determination of Unknown Parameters 


F a (a,t) = a( a)e hJ 


(18) 


where F a (a f t) denotes the deficiency function for the 
unsteady lift or pitching moment coefficient. With this 
assumption, the unsteady lift or pitching moment coef- 
ficient is given by 

y(f) = -j a(a(x))e~ hil ~ x) a(x)dx (19) 

0 

Here, a( a) is an unknown parameter that is a function 
of angle of attack and b i is an unknown constant. Note 
that the parameters u(a) and b, will have different val- 
ues for lift and pitching moment deficiency functions. 
The time derivative of y{t) is given by 


v(0 = ~ 

dt 


( 20 ) 


ja(a(x))^ u T> d(x )dx 

n 

j (~b i )fl(a(T))e' ft|l '' ,, d(T)t/T + a(a(f))a(/) 


or. 


In this study, one set of experimental data was used 
to estimate the unknown parameters in Model I and 
Model II and other sets of available experimental data 
were used to assess the quality of predictions made us- 
ing the estimated parameters. The genetic algorithm 
(GA) was used as a parameter identification algorithm 
because GA provides a convenient environment to es- 
timate the unknown parameters when their functional 
dependence (model structure) is not known a priori. A 
brief description of the GA is given in the following: 

The genetic algorithms are non derivative search 
procedures based on models of processes in natural 
genetics. At each iteration, the search takes place by 
sampling a coding parameter set from a random distri- 
bution, and applying the parameter iterates to a popula- 
tion of function evaluations. The performance of the 
population elements is used to modify the probability 
distribution for sampling in the next iteration. Algo- 
rithms of this class are known to be very robust and not 
particularly vulnerable to termination at local minima. 
More information on GA may be obtained in Refs 1 1 
and 12. 


v(t) = -b i y(t)-a(a(t))a(t) (2!) 

Thus, Eq. (21) is an equivalent ordinary differential 
equation for the convolution integral form (Eq.(I9)) of 
the unsteady lift coefficient y(f). Given the initial value 
y(0), this differential equation is much more convenient 
to solve compared to the evaluation of the convolution 
integral in Eq.( 19). 

If we assume b x to be a function of angle of attack, 
the convolution integral for y(r) which still has the ordi- 
nary differential equation of the form given in Eq.(21) 
is given by the following expression: 9 

j 

y(t) = -J a(a(x))e r a (x)dx (22) 

o 

We note that the convolution integral in Eq.(22) is more 
complex than the convolution integral in Eq.(14) but we 
need not evaluate it. Instead, we can evaluate y(t) more 
easily by solving the ordinary differential equation as 
given in Eq. (21) with b { replaced by b\(a). In this 
study, the functional forms for these two parameters are 
not prescribed a priori. Instead, such functional forms 
are determined during the estimation process. 


The procedure used to determine the unknown pa- 
rameters is as follows: First, subdivide the angle of at- 
tack range of interest into a convenient number of seg- 
ments or nodes, say /t, and use a linear interpolation to 
determine the parameter value in between the nodal 
points. In GA environment, the sub intervals need not 
be equal. The next step is to specify suitable values for 
upper and lower limits at each nodal point for a*( a), 
a 2 *(a), b ,(a), b 2 ( a), M a) and f° r Model I and 

a( a) and b^a) for Model II and upper and lower limits 
for other unknown parameters such as c h c 2 , c 3 , d x , d 2 , 
d y for both models. The GA generates a population 
whose elements are assigned randomly generated val- 
ues for each of the unknown parameters that lie within 
the specified upper and lower limits. Thus, having as- 
signed values to all unknown parameters in Model I or 
Model II, the time histories of C L {t) or C,„(t) are gener- 
ated for each population element. Let the estimated 
time history be denoted by C L (t ) or C ffl (t). Then define, 

J = i\C a (t i )-C u (t l )] (23) 
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where, C a = C L or C w( , and, = L...N correspond to the 
values of time at which the oscillatory data on C L (t ) or 
Cjt) is available. The problem is to estimate the nodal 
values of the unknown parameters in Model I or Model 
IT so that the sum of the squares of the error J is mini- 
mized. The GA was used as a least square minimization 
algorithm with J as the cost function. 

In the longitudinal oscillatory tests of the F-I6XL 
aircraft model, the angle of attack varied from 0 to a 
little over 70 deg. To cover this range of angle of at- 
tack, we use n = 37 (Acx = 2 deg) for Model I and 
n= 16 (Act = 5 deg) for Model II. More number of 
nodes were used for Model I for better definition of the 
parameter variations. Further, we have to specify cer- 
tain GA specific parameters such as population size, 
probability of crossover and probability of mutation. In 
this study, we use population size 30. probability of 
crossover 0.95 and probability of mutation 0.01. If the 
GA reduces cost as the iterations progress, then the GA 
was allowed to continue. Otherwise, the GA was 
stopped, the values of lower and upper limits were rede- 
fined and the GA was started again. Since the GA has 
no defined convergence criterion to end the iteration 
process, the GA was stopped when it became evident 
that further continuation would not result in any im- 
provement in J. 

The convolution integral was evaluated using the 
trapezoidal rule with a step size of 0.0001 s. The inte- 
gration of the ordinary differential equation was done 
using a forth order Runge-Kutta method with a fixed 
step size of 0.02 s. All the computations were done in 
MATLAB. 

Results and Discussion: 

The above formulation for evaluating unsteady lift 
and moment coefficients based on indicial response 
method was applied to the analysis of longitudinal os- 
cillatory data of the 0.1 scale model of the F-16XL air- 
craft obtained in the Langley 12-Foot Low Speed Wind 
Tunnel. A three-view sketch of the 0.1 scale model of 
the F-16-XL aircraft is presented in Fig. 1. Using this 
sketch, the center of the planform area x was found to 
be approximately at 1.0090 m. The moment reference 
point is located at a distance of 0.9932 m from the nose. 
The mean aerodynamic chord (mac) of the 0.1 scale 
model of the F-16XL aircraft is equal to 0.753 m. Ad- 
ditional information on the test model and test facility 
may be found in Refs. 9 and 1 1 . 



Figure 1. Three-View sketch of the F-16XL aircraft 
model. 


In the oscillatory tests, the mean angle of attack 
was close to 35 deg and the oscillation amplitude was 
also about 35 deg so that the angle of attack varied from 
0 to a little more than 70 deg during one oscillation. 
The test section speed was around 17.5 m/sec that cor- 
responds to a Mach number of about 0.05. The period 
of oscillation T f> varied from 12 s to 1.33 s so that the 
corresponding reduced frequency k was in the range 
0.01 1 to 0.101. For the purpose of parameter estima- 
tion, the oscillatory wind tunnel data corresponding to 
T p = 2.38 s (k = 0.056) was used. The rest of the avail- 
able oscillatory test data were used for validation of 
parameter estimates. 

The variation of the lift and pitching moment coef- 
ficients from static wind tunnel tests on the F-16XL 
aircraft model are presented in Fig. 2. Using numerical 
curve fit to these data points, the lift-curve-slope 
and slope of the pitching-moment-coefficient C^were 



Figure 2. Static wind tunnel data of the F-16XL 
aircraft model. 
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deduced as shown in Fig. 3. However, the values of C Ut 
and C„ JU deduced by two different methods of curve 
fitting such as cubic splines and cubic polynomials dif- 
fered by as much ±10% for C ta and about ±15% for 
C ma , especially at high angles of attack. In view of this, 
we assume 0.9 < a*( a) < 1.1 in Eq. (11) and 0.85 < 
a*((x) < 1.15 in Eq. (12). For better definition of these 
error bounds, more data points taken at much smaller 
intervals are needed. However, such data were not 
available for this study. 



Figure 3. Lift-Curve-Slope and slope of pitching- 
moment-coefficient of the F-16XL aircraft model. 

The unknown coefficients in Eq. (15) - (17) for C U/ 
and C mq were estimated as follows: Model I: 
Cj = -0.3844, c 2 = 2.3776, c, = -4.5026, d l = -0.3678, 
d 2 = -0.1616, and d 3 = 2.0361; For Model II: 
c] = -0.4240, c 2 = 3.3127, c 3 = -3.3840, d, = -1.2450, 
d 2 = -0.3806, and d 3 = 1 .5557. 



0 2 4 6 8 10 12 

Time, sec 


Figure 4. Time history of angle of attack, T p = 2.38 s. 

1.6 
1.4 
1.2 
1.0 
C L 0.8 
0.6 
0.4 
0.2 

0 2 4 6 8 10 12 

Time, sec 

Figure 5. Time history 7 of measured and estimated 
lift coefficient, T p = 2.38 s. 


° Oscillatory WT Data 
- - - Mode! ! 

Model II 



The lime history of the angle of attack tor 
T p = 2.38 s (k = 0.056) is shown in Fig. 4. The esti- 
mated and measured time histories of lift coefficient are 
shown in Fig. 5. A cross plot of C L (t) with a(/) is shown 
in Fig. 6. It is observed the estimated time histories 
match the measured history fairly well, with Model II 
giving better fit to the oscillatory test data. A cross plot 
of the corresponding time varying lift coefficient with 
time varying angle of attack is shown in Fig. 6. It was 
observed that the contribution of rotary derivative is 
very small and negligible. On the other hand, the con- 
tribution of unsteady term y(t) happens to be signifi- 
cant. It is quite possible that if more accurate values of 
C La were available, Model I may have given better re- 
sult. The estimated values of the indicial parameters 
Model I and Model II are shown in Fig. 7(a) and 7(b). 
For comparison, the lift-curve-slope is included in Fig. 
7(b). It is observed that the estimated values of a 2 *( a) 


a Oscillatory WT Data 

Model I 

Model II 



Figure 6. Comparison of measured and estimated 
lift coefficient, T p - 2.38 s. 
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are in the range 0.38 to 0.46 suggesting that the initial 
peak value of the noncirculatory loading is only 0.38 to 
0.46 of that based on linear piston theory. 



ot, deg 


Figure 7(a) Estimated indicial function parameters 
in lift coefficient, Model I. 



0 10 20 30 40 50 60 70 


a, deg 

Figure 7(b). Estimated indicial function parameters 
in lift coefficient, Model II. 

Using these parameter estimates as shown in Fig. 
7, variations of the noncircular and circular components 
as well total loading for Model I (Eq. (11)), and total 
loading for Model II (Eq.(!8)) were calculated for 
0 < a < 70 deg. Some sample results of these variations 
are shown in Figs. 8 to 10. As expected, the noncircu- 
latory component decays rapidly and the circulatory 
component decays gradually. The area under each of 
these curves is a measure of their contribution to the 
unsteady lift coefficient. The calculated areas are shown 
in Fig. 11. To first order in frequency, the area under 
the deficiency function is related to the parameter 
or C ma as the case may be. Thus, the process of esti- 
mating C Lli and area under the deficiency function is 
equivalent to estimating the combined parameter 



Time, s 


Figure 8. Estimated deficiency function components 
in lift coefficient, T p - 2.38 s, a = 15 deg 



Time, s 


Figure 9. Estimated deficiency function components 
in lift coefficient, T p = 2.38 s, a = 30 deg. 



Time, s 

Figure 10. Estimated deficiency function 
components in lift coefficient, T p - 2.38 s, a = 50 deg. 
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C^ + C^ or C„ tt/ + C /na in small amplitude oscillatory 
wind tunnel tests. From Fig. 1 1 we observe that the 
contribution of noncircular component is always nega- 
tive and gradually drops with increase in angle of at- 
tack. The sign of the circulatory component follows the 
sign of C L (a). It is interesting to note that Model II fol- 
lows closely the variation of the total loading of Model 
I which is the sum of the noncirculatory and circulatory 
components. Thus, the Model II with one term expo- 
nential function satisfactorily approximates the com- 
bined variation of the noncirculatory and circulatory 
components of Model I. Further, Model II has simple 
structure and is computationally efficient because the 
convolution integral in the deficiency function can be 
represented in the form of an ordinary differential 
equation that can be solved more easily compared to the 
evaluation of convolution integrals in Model I. Further, 
the availability of ordinary' differential equation greatly 
facilitates the inclusion of unsteady aerodynamic terms 
in flight simulation computer codes. 



Figure 11. Variation of area under estimated 
deficiency function components, T p = 2.38 s. 

The estimated time histories of the pitching mo- 
ment coefficient T p = 2.38s (k = 0.0562) are shown in 
Fig 12. The corresponding cross plots of the time vary- 
ing pitching moment coefficient with time varying an- 
gle of attack are shown in Figure 13. The estimated 
values of the indicial parameters for Model I and Model 
II are shown in Figures 14(a) and 14(b). We observe 
that the quality of the fit between the oscillatory pitch- 
ing moment coefficient data and the estimates is not as 
good as that obtained for the lift coefficient. In this case 
also, Model II gives a better Fit to the oscillatory data. 

Using these parameter estimates, predictions were 
made for T p = 1.33s, 1. 72s, 4.0s and 12.0s for which 
oscillatory data was obtained. These results are shown 


O Oscillatory WT Data 
— Mode!! 

--- Modelll 



Figure 12. Time history of measured and estimated 
pitching moment coefficient, T p = 2.38 s. 
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Figure 13. Comparison of measured and estimated 
pitching moment coefficient, T p = 2.38 s. 



Figure 14(a). Estimated indicial function parameters 
in pitching moment coefficient, Model I. 
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a, deg 


Figure 14(b). Estimated indicia! function parameters 
in pitching moment coefficient, Model IT. 

in Figs. 15 to 22. It is observed that the predicted time 
histories of lift and moment coefficients agree reasona- 
bly well with measured oscillatory wind tunnel data, 
with better agreement indicated for the lift coefficient. 
In all these cases, Model II gives better predictions 
compared to Model I. However, differences do exist, 
particularly in the pitching moment coefficient. The 
quality of predictions of Model I might have been better 
if more static wind tunnel data at smaller intervals were 
available for better definition of C La > C ma and the upper 
and lower bounds for the uncertainty parameter a y \o). 

The indicial response method of Tobak 2 3 is sup- 
posed to be applicable when the lift-curve-slope is 


° Oscillatory WT Data 

Model ! 

Model II 



Time, sec 


Figure 15. Time history of measured and predicted 
lift coefficient, T p = 1.33 s. 


positive and the steady stale of the indicial response is 
time invariant. However, in the present study, the lifl- 
curve-slopc goes negative beyond 30 deg angle of at- 
tack. Further, it is quite possible that the flow field over 
the F-16XL aircraft model at high angles of attack is 
characterized by extensive flow separations, vortex 
interactions and possible vortex bursts over wing and 
tail surfaces. Under such flow conditions, the assump- 
tion of time invariant steady stales may not be strictly 
valid and the analysis should include the possible exis- 
tence of time-dependent equilibrium states as suggested 
by Tobak et al, 54 which was not attempted in this study. 


o Oscillatory WT Data 
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Model II 
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Time, sec 

Figure 16. Time history of measured and predicted 
lift coefficient, T p = 1.72 s. 
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Figure 17. Time history of measured and predicted 
lift coefficient, T p = 4.0 s. 
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o Oscillatory WT Data 
Model 1 



Time, sec 

Figure 18. Time history of measured and predicted 
lift coefficient, T p = 12.0 s. 
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Figure 19. Time history of measured and predicted 
pitching moment coefficient, T p = 1.33 s. 



Figure 20. Time history of measured and predicted 
pitching moment coefficient, T p = 1.72 s. 
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Figure 21. Time history of measured and predicted 
pitching moment coefficient, T p = 4.0 s. 
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Figure 22. Time history of measured and predicted 
pitching moment coefficient, T p = 12.0 s. 

Concluding Remarks 

The nonlinear indicial response method provides a 
fundamental approach to the problem of modeling un- 
steady aerodynamic effects at high angles of attack. In 
this study, two models have been used to approximate 
the nonlinear indicial response. Model I is based on two 
term exponential functions, with one term to represent 
the noncirculatory component and the other term to 
represent the decay of circulatory component. Some 
adjustments have been introduced to Model I to account 
for uncertainties in the static wind tunnel data that is 
needed in Model I and for deviations from linear piston 
theory due to nonlinearities at high angles of attack. 
Model II is based on one exponential function to ap- 
proximate the combined variation of noncirculatory and 
circulatory components. The unknown indicial pa- 
rameters in each of the two models were estimated us- 
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ing a genetic algorithm. The genetic algorithm was used 
as a least-square minimizing technique. 

It is shown that Model II is capable of satisfactorily 
approximating the combined variation of the noncircu- 
latory and circulatory components of Model I. Further, 
Model II is simple in structure and computationally 
efficient because it can be expressed in the form of an 
ordinary differential equation as opposed to the convo- 
lution integrals in Model I. This fact greatly facilitates 
the inclusion of unsteady aerodynamic terms in the 
flight simulation computer codes. 
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